
/*
* Additional data prep for Tab5 and Fig5, using Indice de Marginacion
*/

*Pollution from fires
use "$dataDir/dataset_pollution_harvest_daily_weather", clear
		
		local depvar "pm25"
		local cluslvl "id_mill"
		
foreach var in `depvar'{
	reghdfe `var' afterXtreatment temp_celsius rain rain2, abs(id_mills i.year#i.month) vce(cluster `cluslvl')
		predict `var'_pred, xb
	}
	
	rename date date_dmy
	collapse (mean) pm25_pred temp_celsius, by(id_mill date_dmy after2015 id_mun treatment year month)

	keep id_mills id_mun year month date_dmy pm25_pred after2015
	gen id=id_mills
	
	merge 1:m id date_dmy using "$dataDir/dataset_dac_locs_catchmentareas_fires", gen(mergedac)
		
*DISADVANTAGED COMMUNITIES DEFINITION	
	gen DAC=(gm_2010=="Alto" | gm_2010=="Muy alto")
*GROUPS
	gen group=0 if gm_2010=="Muy bajo"
		replace group=1 if gm_2010=="Bajo"
		replace group=2 if gm_2010=="Medio"
		replace group=3 if gm_2010=="Alto"
		replace group=4 if gm_2010=="Muy alto"
	
	encode id_loc, gen(idloc)
	
	gen year2=year(date_dmy)
	gen month2=month(date_dmy)


	collapse (mean) pm25_pred pob_tot, by(id_mill after2015 id_mun after2015 year DAC group idloc)
	reghdfe pm25_pred i.DAC##after2015, abs(year idloc) vce(cluster idloc)
		estimates store model1
		qui sum pm25_pred if e(sample)==1 & year<2015
		local mean1=r(mean)
		estadd scalar mean1 `mean1' : model1 
		estadd local yearfe "Yes"
		estadd local locfe "Yes"
		estadd local cluslvl "Locality"
	reghdfe pm25_pred i.group##after2015, abs(year idloc) vce(cluster idloc)
		estimates store model1_levels
		
*Pollution from mills	

use "$dataDir/dataset_pollution_mills_year_weather", replace

		local depvar "pm25"
		local cluslvl "id_mun year"
		
foreach var in `depvar'{

	reghdfe `var' afterXtreatment temp_celsius rain rain2, abs(id_mills i.year) vce(cluster `cluslvl')

			predict `var'_pred, xb
	}
	
	
	joinby id_mills using "$dataDir/dataset_dac_locs_catchmentareas_mills", unm(b) 
	
*DISADVANTAGED COMMUNITIES DEFINITION	
	gen DAC=(gm_2010=="Alto" | gm_2010=="Muy alto")
*GROUPS
	gen group=0 if gm_2010=="Muy bajo"
		replace group=1 if gm_2010=="Bajo"
		replace group=2 if gm_2010=="Medio"
		replace group=3 if gm_2010=="Alto"
		replace group=4 if gm_2010=="Muy alto"
		
	encode id_loc, gen(idloc)
	

	reghdfe pm25_pred i.DAC##after2015, abs(year idloc) vce(cluster idloc)
		estimates store model2
		qui sum pm25_pred if e(sample)==1 & year<2015
		local mean1=r(mean)
		estadd scalar mean1 `mean1' : model2 
		estadd local yearfe "Yes"
		estadd local locfe "Yes"
		estadd local cluslvl "Locality"
			
	reghdfe pm25_pred i.group##after2015, abs(year idloc) vce(cluster idloc)
		estimates store model2_levels

	coefplot model1_levels model2_levels, omitted baselevels keep(0.group#1.after2015 1.group#1.after2015 2.group#1.after2015 3.group#1.after2015 4.group#1.after2015) coeflabel(0.group#1.after2015 = "Very low" 1.group#1.after2015 = "Low" 2.group#1.after2015 = "Medium" 3.group#1.after2015 = "High" 4.group#1.after2015 = "Very High") order(0.group#1.after2015 1.group#1.after2015 2.group#1.after2015 3.group#1.after2015 4.group#1.after2015) vertical yline(0) xtitle("Disadvantaged community index") ciopts(recast(rcap)) ytitle("Change on `poll' ({&mu}g/m{sup:3}) concentration after 2015", size(small)) p1(label("Pollution from fires"))  p2(label("Pollution from mills"))
			graph export "$mainDir/figures/Fig5_panelA_distributional_outcomes_pm25.pdf", as(pdf) replace		

			esttab model1 model2 using "$mainDir/tables/Tab5_distributional_outcomes_pm25.tex", replace b(%12.3fc) se(%12.3fc) keep(1.DAC#1.after2015) varlabels(1.after2015#1.treatment  "After 2015 $\times$ Disadvantaged") label starl(* 0.1 ** 0.05 *** 0.01) stats(mean1 N r2 yearfe locfe cluslvl, fmt(%9.3fc %12.0fc %12.3fc) label("Pre 2015 mean" "Obs." "R-squared" "Year FE" "Locality FE" "Cluster level")) nonotes		

************************************************************************************************************************************			
/*
* Additional data prep for Tab5 and Fig5, using Indice de Rezago Social
*/

*Pollution from fires
use "$dataDir/dataset_pollution_harvest_daily_weather", clear
		
		local depvar "pm25"
		local cluslvl "id_mill"
		
foreach var in `depvar'{
	reghdfe `var' afterXtreatment temp_celsius rain rain2, abs(id_mills i.year#i.month) vce(cluster `cluslvl')
			predict `var'_pred, xb
	}
	
	rename date date_dmy
	
	collapse (mean) pm25_pred temp_celsius, by(id_mill date_dmy after2015 id_mun treatment year month)
	keep id_mills id_mun year month date_dmy pm25_pred after2015
		gen id=id_mills

	merge 1:m id date_dmy using "$dataDir/dataset_rezago_locs_catchmentareas_fires", gen(mergedac)
	
		
	
*DISADVANTAGED COMMUNITIES DEFINITION	
	gen DAC=(grado_rezago_social=="Alto" | grado_rezago_social=="Muy alto")
	*GROUPS
	gen group=0 if grado_rezago_social=="Muy bajo"
		replace group=1 if grado_rezago_social=="Bajo"
		replace group=2 if grado_rezago_social=="Medio"
		replace group=3 if grado_rezago_social=="Alto"
		replace group=4 if grado_rezago_social=="Muy alto"
	
	encode id_loc, gen(idloc)
	
	gen year2=year(date_dmy)
	gen month2=month(date_dmy)
	
	rename total_pop pob_tot
	collapse (mean) pm25_pred pob_tot, by(id_mill after2015 id_mun after2015 year DAC group idloc)
	reghdfe pm25_pred i.group##after2015, abs(year idloc) vce(cluster idloc)
		estimates store model1_levelsr
		
*Pollution from mills	

use "$dataDir/dataset_pollution_mills_year_weather", replace
			
		local depvar "pm25"
		local cluslvl "id_mun year"

foreach var in `depvar'{
	reghdfe `var' afterXtreatment temp_celsius rain rain2, abs(id_mills i.year) vce(cluster `cluslvl')
			predict `var'_pred, xb
	}
	
	
	joinby id_mills using "$dataDir/dataset_rezago_locs_catchmentareas_mills"
	
*DISADVANTAGED COMMUNITIES DEFINITION	
	gen DAC=(grado_rezago_social=="Alto" | grado_rezago_social=="Muy alto")
	*GROUPS
	gen group=0 if grado_rezago_social=="Muy bajo"
		replace group=1 if grado_rezago_social=="Bajo"
		replace group=2 if grado_rezago_social=="Medio"
		replace group=3 if grado_rezago_social=="Alto"
		replace group=4 if grado_rezago_social=="Muy alto"
	
	encode id_loc, gen(idloc)
		rename total_pop pob_tot
		
reghdfe pm25_pred i.group##after2015, abs( year idloc) vce(cluster idloc)
	estimates store model2_levelsr

***************************************************************************************************************

	coefplot model1_levelsr model2_levelsr, omitted baselevels keep(0.group#1.after2015 1.group#1.after2015 2.group#1.after2015 3.group#1.after2015 4.group#1.after2015) coeflabel(0.group#1.after2015 = "Very low" 1.group#1.after2015 = "Low" 2.group#1.after2015 = "Medium" 3.group#1.after2015 = "High" 4.group#1.after2015 = "Very High") order(0.group#1.after2015 1.group#1.after2015 2.group#1.after2015 3.group#1.after2015 4.group#1.after2015) vertical yline(0) xtitle("Disadvantaged community index") ciopts(recast(rcap)) ytitle("Change on `poll' ({&mu}g/m{sup:3}) concentration after 2015", size(small)) p1(label("Pollution from fires"))  p2(label("Pollution from mills"))
			graph export "$mainDir/figures/Fig5_panelB_distributional_outcomes_pm25.pdf", as(pdf) replace		
